Anaerobic CUE

set.seed(1000)
samples_prep_path <- "./"
data_path <- "./DADA2_pseudo/"
Metadata_table <- "./AnCUE_Metadata.csv"
Seq_table <- "DADA2.seqtab_nochim.tsv"
Tax_table <- "DADA2.taxa_silva.tsv"
Seq_file <- "DADA2.Seqs.fa"

# read OTU mat from data file
read_tsv(paste0(data_path, Seq_table), 
                        trim_ws = TRUE) %>% 
  column_to_rownames("ASV") %>% 
  t() %>% 
  as.data.frame() -> # not tibble because we need row names
  abundance_mat # convert to abundance matrix

# get short names of samples
# abundance_mat %>% 
#   rownames() %>% 
#   str_remove("^Roey[0-9]{3,}-?") %>% 
#   str_split("_", simplify = T) %>% 
#   .[, 1] ->
#   short_names

# Read metadata file
read_csv(paste0(samples_prep_path, Metadata_table), 
         trim_ws = TRUE) %>% 
  mutate(`16S copies` = replace(`16S copies`, which(`16S copies` == 0 | is.na(`16S copies`)), 1)) %>%  # add pseudo count
  filter(merged_sample_name %in% str_remove(rownames(abundance_mat), "_L001")) %>% # remove metadata rows if the samples did not go through qual processing
  mutate(`Library size` = rowSums(abundance_mat)) %>% # Add lib size
  mutate(to_names = merged_sample_name) %>% 
  mutate(across(c(
    "Site",
    "Oxygen",
    "Glucose",
    "Label (13C)"), 
    ~factor(.))) %>% 
  mutate(`Density zone` = factor(ifelse(`Density (g ml-1)` > 1.795, "Heavy", "Light"), levels = c("Light", "Heavy"))) %>% # critical for DESeq2 that the reference is the first level
  column_to_rownames("to_names") ->
  Metadata

# Order abundance_mat samples according to the metadata
sample_order <- match(str_remove(rownames(abundance_mat), "_L001"), rownames(Metadata))
abundance_mat %<>% arrange(sample_order)
rownames(abundance_mat) <- rownames(Metadata) # needed for pyhloseq

# read taxonomy from data file
Raw_tax_data <- read_tsv(paste0(data_path, Tax_table), 
                        trim_ws = TRUE, col_names = TRUE)
Raw_tax_data %<>%
  mutate_all(~(replace(., is.na(.), "Unclassified"))) # I think mutaute_all is unnecessary here because replace(., is.na(.), "Unclassified") alone should work

Raw_tax_data %>%
  dplyr::select(.,
         `Kingdom (BS)`,
         `Phylum (BS)`,
         `Class (BS)`,
         `Order (BS)`,
         `Family (BS)`,
         `Genus (BS)`) %>%
  cbind(Name = colnames(abundance_mat),. ) ->
  Taxonomy.bs

Raw_tax_data %>%
  dplyr::select(.,
         Kingdom,
         Phylum,
         Class,
         Order,
         Family,
         Genus) %>% 
  # map_dfr(., as_factor) %>% 
  # map_dfr(fct_expand, "Rare")  %>%  
  as.matrix() -> # must be a matrix or phyloseq drops row names and gives and error
  Taxonomy
row.names(Taxonomy) <- colnames(abundance_mat)
# colnames(Taxonomy) <-
#   c("Domain", "Phylum", "Class", "Order", "Family", "Genus")

# read sequence data
ASV_seqs <- readDNAStringSet(
  file = paste0(data_path, Seq_file),
  format = "fasta", 
  nrec = -1L, 
  skip = 0L, 
  seek.first.rec = FALSE,
  use.names = TRUE)

# generate phyloseq object
Ps_obj <- phyloseq(otu_table(abundance_mat, taxa_are_rows = FALSE),
                   sample_data(Metadata),
                   tax_table(Taxonomy),
                   refseq(ASV_seqs))

Ps_obj %>% 
  get_variable() %>% 
  vis_dat()

Ps_obj %>% 
  get_variable() %>% 
  select_if(is.numeric) %>% 
  vis_cor()

Ps_obj_df <-
  as.data.frame(sample_data(Ps_obj)) # Put sample_data into a ggplot-friendly data.frame
Ps_obj_df <- Ps_obj_df[order(Ps_obj_df$`Library size`), ]
Ps_obj_df$Index <- seq(nrow(Ps_obj_df))
ggplot(data = Ps_obj_df, 
       aes(x = Index, y = Library.size, color = Control)) + 
  geom_point() +
  scale_y_log10(breaks = c(
    min(Ps_obj_df$Library.size),
    10,
    100,
    1000,
    5000,
    10000,
    ceiling(max(Ps_obj_df$Library.size) / 10000) * 10000
    )) + 
  scale_color_brewer(type = 'qual', palette = 'Set1', direction = -1)

summary(sample_sums(Ps_obj))

summary(taxa_sums(Ps_obj))

Ps_obj %<>%
  prune_samples(names(which(sample_sums(Ps_obj) != 0)), .)
# Ps_obj <-
#   subset_samples(Ps_obj, Sample != "CTRL")
# summary(sample_sums(Ps_obj))

contamdf.freq <-
  isContaminant(Ps_obj, method = "frequency", conc = "X16S.copies")
# print(contamdf.freq)
# How many contaminants are found?
table(contamdf.freq$contaminant)

# Which ones
which(contamdf.freq$contaminant)

plot_frequency(Ps_obj, taxa_names(Ps_obj)[c(1, 3)], conc = "X16S.copies")

plot_frequency(Ps_obj, taxa_names(Ps_obj)[which(contamdf.freq$contaminant)[1:20]], conc = "X16S.copies")

contamdf.prev <- isContaminant(Ps_obj, method = "prevalence", neg = "Control")
# How many contaminants are found?
table(contamdf.prev$contaminant)

# Which ones
which(contamdf.prev$contaminant)

# And using a more aggressive threshold
contamdf.prev05 <- isContaminant(Ps_obj, method = "prevalence", neg = "Control", threshold = 0.5)
table(contamdf.prev05$contaminant)

# Make phyloseq object of presence-absence in negative controls
Ps_obj.pa <-
  transform_sample_counts(Ps_obj, function(abund)
    1 * (abund > 0))
Ps_obj.pa.neg <-
  prune_samples(sample_data(Ps_obj.pa)$Control == "TRUE", Ps_obj.pa)
Ps_obj.pa.pos <-
  prune_samples(sample_data(Ps_obj.pa)$Control == "FALSE", Ps_obj.pa)
# Make data.frame of prevalence in positive and negative samples
df.pa <-
  data.frame(
    pa.pos = taxa_sums(Ps_obj.pa.pos),
    pa.neg = taxa_sums(Ps_obj.pa.neg),
    contaminant = contamdf.prev$contaminant
  )
ggplot(data = df.pa, aes(x = pa.neg, y = pa.pos, color = contaminant)) + geom_point() +
  xlab("Prevalence (Negative Controls)") + ylab("Prevalence (True Samples)")

c(taxa_names(Ps_obj)[which(contamdf.freq$contaminant)],
  taxa_names(Ps_obj)[which(contamdf.prev$contaminant)]) ->
  contaminant_seqs
  
write_csv(as_tibble(contaminant_seqs), 
            paste0(data_path, "decontam_contaminants.csv"), 
            col_names = FALSE)


good_seqs <- setdiff(taxa_names(Ps_obj), contaminant_seqs)
Ps_obj_clean <- prune_taxa(good_seqs, Ps_obj)

# save decontaminated seqtab
Ps_obj_clean %>% 
  t() %>% 
  get_taxa() %>%
  as_tibble(rownames = "ASV") %>%
  write_tsv(., 
            paste0(data_path, str_remove(Seq_table, ".tsv"), "_decontam.tsv"), 
            col_names = TRUE)

Ps_obj_clean %>% 
  t() %>% 
  tax_table() %>%
  as_tibble(rownames = "ASV") %>%
  write_tsv(., 
            paste0(data_path, str_remove(Tax_table, ".tsv"), "_decontam.tsv"), 
            col_names = TRUE)

# save decontaminated metadata (just in case some samples were dropped)
Ps_obj_clean %>% 
  t() %>% 
  get_variable() %>% 
  setNames(., colnames(Metadata)) %>% 
  # as_tibble(rownames = "ASV") %>%
  write_csv(., 
            paste0("./", str_remove(Metadata_table, ".csv"), "_decontam.csv"), 
            col_names = TRUE)

# save decontaminated seqs
Ps_obj_clean %>% 
   refseq() %>% 
  writeXStringSet(., filepath = paste0(data_path, str_remove(Seq_file, ".fa*"), "_decontam.fa"), format = "fasta", width = 1000)
 
# save R obj
saveRDS(Ps_obj_clean, file = paste0(data_path, "Ps_obj_decontam.Rds"))

sessioninfo::session_info() %>%
  details::details(
    summary = 'Current session info',
    open    = TRUE
  )